#include <RcppArmadillo.h>
#include <RcppArmadilloExtensions/sample.h>

// [[Rcpp::depends(RcppArmadillo)]]
# include <omp.h>

// [[Rcpp::plugins(openmp)]]

// state vars
// > data.frame(pos = 0:(length(state.vars)-1), var = state.vars)
//   pos                  var
// 1   0                htype
// 2   1                 size
// 3   2                 rent
// 4   3          stay_length
// 5   4          gwealth_imp
// 6   5                  age
// 7   6               income
// 8   7           two_adults
// 9   8                 kids
// 10  9         rent_htype01
// ...
// 19 18         rent_htype10


// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::plugins("cpp11")]]


// compute income transitions
void inc_transition(
    Rcpp::NumericVector& inc_lead,
    Rcpp::NumericVector& inc_transition_coefs,
    Rcpp::NumericVector& inc_transition_errors,
    Rcpp::NumericVector& state,
    const int& two_adults_lead
) {

  int K = inc_transition_coefs.length();
  Rcpp::NumericVector vars(K);

  vars(0) = 1; // Intercept
  vars(1) = std::log(state(6)); // log income
  vars(2) = state(5); // age
  vars(3) = std::pow(state(5), 2) * 1e-3; // age squared
  vars(4) = state(7); // two_adults
  vars(5) = two_adults_lead; // two_adults_lead

  double V = 0 ;
  for(int k = 0; k < K; k++) {
    V += vars(k) * inc_transition_coefs(k);
  }

  int E = inc_transition_errors.length();
  for(int e = 0; e < E; e++) {
    inc_lead(e) = std::exp(V + inc_transition_errors(e));
  }

  return ;
}

// Compute wealth transitions
void sav_transition(
    Rcpp::NumericVector& sav_lead,
    Rcpp::NumericVector& sav_transition_coefs,
    Rcpp::NumericVector& sav_transition_errors,
    Rcpp::NumericVector& state,
    Rcpp::NumericVector& inc_lead,
    const int& j0,
    Rcpp::NumericVector& sizes_sqm
) {
  int K = sav_transition_coefs.length();
  Rcpp::NumericVector vars(K);

  // vars(0) = 1 ; // intercept
  // vars(1) = 0 ; // expected income
  vars(2) = (state(2) + 2.5*sizes_sqm(state(1) - 1))*12; // rent + other housing costs
  vars(3) = (j0 < 0) ? 0 : 1;  // move indicator
  vars(4) = state(4) ; // gwealth_imp

  double V = sav_transition_coefs(0) ; // intercept
  for(int k = 2; k < K; k++) {
    V += vars(k) * sav_transition_coefs(k);
  }
  int E = sav_transition_errors.length();
  for(int k = 0; k < E; k++) { // savings depend on realized future income + the error
    sav_lead(k)  = V + inc_lead(k) * 12 * sav_transition_coefs(1); // expected income, t+1
    sav_lead(k) += sav_transition_errors(k);
    sav_lead(k)  = (sav_lead(k) < 0) ? 0 : sav_lead(k); // no net debt allowed
  }
  return ;
}

//
// couple indicator transition
void adu_transition_prob(
    Rcpp::NumericVector& prob,
    Rcpp::NumericMatrix& adu_transition_coefs,
    Rcpp::NumericVector& state
) {
  int col = state(7); // two adults (select model in col = 0, 1)
  int K = adu_transition_coefs.nrow();
  Rcpp::NumericVector vars(K);

  vars(0) = 1;
  vars(1) = state(5); // age
  vars(2) = std::pow(state(5), 2) * 1e-3; // age squared
  vars(3) = (state(8) == 1) ? 1 : 0; // 1 kid
  vars(4) = (state(8) == 2) ? 1 : 0; // 2+ kids

  double V = 0;
  for(int k = 0; k < K; k++) {
    V += adu_transition_coefs(k,col) * vars(k);

  }
  V = std::exp(V);

  prob(1) = V / (1 + V) ;
  prob(0) = 1 - prob(1) ;

  return ;
}

//
// kids transitions
void kid_transition_prob(
    Rcpp::NumericVector& prob,
    Rcpp::NumericMatrix& kid_transition_coefs,
    Rcpp::NumericVector& state
) {

  int K = kid_transition_coefs.nrow();
  Rcpp::NumericVector vars(K);

  vars(0) = 1;
  vars(1) = state(5); // age
  vars(2) = std::pow(state(5), 2) * 1e-3; // age squared
  vars(3) = state(7); // two adults
  vars(4) = (state(8) == 1) ? 1 : 0; // 1 kid
  vars(5) = (state(8) == 2) ? 1 : 0; // 2+ kids

  Rcpp::NumericVector V(2);
  double V_sum = 1;
  for(int i = 0; i < 2; i++) {
    V(i) = 0;
    for(int k = 0; k < K; k++) {
      V(i) += kid_transition_coefs(k,i) * vars(k);
    }
    V(i) = std::exp(V(i));
    V_sum += V(i);
  }

  prob(0) = 1 / V_sum ;
  prob(1) = V(0) / V_sum ;
  prob(2) = V(1) / V_sum ;

  return ;
}



// Update determinstic part of state
void update_det_state(
    Rcpp::NumericVector& new_state,
    Rcpp::NumericVector& state,
    const int& j0,
    Rcpp::NumericVector& sizes_sqm,
    const double& g
) {

  for(int i = 0; i < 9; i++) { // copy state
    new_state(i) = state(i);
  }
  for(int i = 9; i < state.length(); i++) { // copy + update rent levels
    new_state(i) = state(i) * (1 + g);
  }

  new_state(5)++; // age + 1

  if(j0 < 0) { // stay
    new_state(3)++; // stay length
  } else { // move
    new_state(1) = j0 / 10 + 1; // size
    new_state(0) = j0 - (new_state(1) - 1)*10 + 1; // htype
    new_state(3) = 0; // stay length
    new_state(2) = sizes_sqm(new_state(1) - 1) * new_state(8 + new_state(0)) ; // rent
  }

  return ;
}

// Same, but for simulation (allowing for market-segment specific rents)
void update_det_state_hs(
    Rcpp::NumericVector& new_state,
    Rcpp::NumericVector& state,
    const int& j0,
    Rcpp::NumericVector& sizes_sqm,
    Rcpp::NumericMatrix& rent_hs_new,
    const double& g
) {
  const int size = j0 / 10 + 1;


  for(int i = 0; i < 9; i++) { // copy state
    new_state(i) = state(i);
  }
  for(int i=0; i<rent_hs_new.ncol(); i++) {
    new_state(9 + i) = rent_hs_new(size - 1, i) * (1 + g); // rent/sqm by htype
  }

  new_state(5)++; // age + 1

  if(j0 < 0) { // stay
    new_state(3)++; // stay length
  } else { // move
    new_state(1) = size; // size
    new_state(0) = j0 - (new_state(1) - 1)*10 + 1; // htype
    new_state(3) = 0; // stay length
    new_state(2) = sizes_sqm(new_state(1) - 1) * new_state(8 + new_state(0)) ; // rent
  }

  return ;
}

// Overall utility function
void util(
    Rcpp::NumericVector& mapped, // storage for utility values
    Rcpp::NumericVector& state, // current state
    const int& j0, // choice
    const int& S, // number of unobserved states
    const double& sw, // signed weight
    const int& nvars_u_c, // number of non-het utility variables
    const int& nvars_u_s, // number of het utility variables
    const int& nvars_v, // number of terminal utility variables
    Rcpp::NumericVector& sizes_sqm // size/sqm for units w/ 1, 2, 3, 4 rooms
) {
  int j = j0 + 1;

  double two_adults = state(7);
  double kids = state(8);

  if(j == 0) {
    int size = state(1);
    // int size_diff = size - (1 + two_adults +  kids) ;
    double log_size_ratio = std::log(size) - std::log(1 + two_adults +  kids) ;

    double disp_inc = state(6) - state(2) - 2.5*sizes_sqm(size - 1) ; // income - rent
    // 2.5 EUR/sqm costs for utilities and heating
    //  Based on: "Pressemitteilung: Neuer Betriebskostenspiegel für Deutschland", from Nov. 30, 2011.
    //  Values refer to 2009.
    //  https://web.archive.org/web/20130726165023/http://www.mieterbund.de/pressemitteilung.html?&no_cache=1&tx_ttnews%5Btt_news%5D=326&tx_ttnews%5BbackPid%5D=3126&cHash=11b626ca54
    // Value for 2009: 2.28 EUR/sqm (2009 EUR)
    // deflated to 2018 EUR: 2.19 / 87.2 * 98.1 =~ 2.5
    //


    disp_inc = disp_inc / sqrt(1 + two_adults +  0.5*kids) ; // equivalized (dep. on hh size)
    double hqual = state(0) * std::exp(-0.038 * state(3));

    // Utility function values
    mapped(0) += (disp_inc < .1) ?  -2.302585 * sw : std::log(disp_inc) * sw ; // income

    mapped(1) += hqual * sw; // hqual

    mapped(2) += sw * log_size_ratio;
    mapped(3) += state(3) * sw ; // stay length

  } else {
    // moving costs
    mapped(nvars_u_c) += sw; // move dummy
    mapped(nvars_u_c + 1) += state(5) * sw; // variable moving costs (age)
    mapped(nvars_u_c + 2) += std::pow(state(5),2) * sw; // variable moving costs (age^2)
    mapped(nvars_u_c + 3) += two_adults * sw; // two adults
    mapped(nvars_u_c + 4) += kids * sw; // number of kids

    if(j <= 40) { // resets stay length to zero + rent & size & htype adjust

      int size = j0 / 10 + 1; // size
      // int size_diff = size - (1 + two_adults +  kids) ;
      double log_size_ratio = std::log(size) - std::log(1 + two_adults +  kids) ;

      int hqual = j - (size - 1)*10; // htype
      double rent = sizes_sqm(size - 1) * state(8 + hqual);
      // double rent = sizes_sqm(size - 1) * rent_hs(size - 1, 8 + hqual); // rent

      double disp_inc = state(6) - rent - 2.5*sizes_sqm(size - 1); // income - rent
      disp_inc = disp_inc / sqrt(1 + two_adults +  0.5*kids) ; // equivalized (dep. on hh size)

      // Utility function values
      mapped(0) += (disp_inc < .1) ?  -2.302585 * sw : std::log(disp_inc) * sw ; // income

      mapped(1) += hqual * sw; // hqual

      mapped(2) += sw * log_size_ratio ;
      // mapped(3) += state(3) * sw // stay length = 0

    } else {

      int select = nvars_u_c + nvars_u_s + (j - 41) * nvars_v ;
      mapped(select + 0) += sw; // constant
      mapped(select + 1) += std::log(state(6)) * sw; // log income
      mapped(select + 2) += state(4)*1e-5 * sw; // wealth
      mapped(select + 3) += std::pow(state(4)*1e-5,2) * sw; // wealth^2
      mapped(select + 4) += state(5) * sw ; // age
      mapped(select + 5) += std::pow(state(5),2) * sw ; // age^2
      mapped(select + 6) += two_adults * sw; // two adults
      mapped(select + 7) += kids * sw; // number of kids

    }
  }

  return ;
}

// utility function used in simulation
void util_hs( // see function util() for comments
    Rcpp::NumericVector& mapped,
    Rcpp::NumericVector& state, // current state
    const int& j0, // current choice
    const int& S, // number of unobserved states
    const double& sw, // signed weight
    const int& nvars_u_c,
    const int& nvars_u_s,
    const int& nvars_v,
    Rcpp::NumericVector& sizes_sqm,
    Rcpp::NumericMatrix& rent_hs // rents in each segment
) {
  int j = j0 + 1;

  double two_adults = state(7);
  double kids = state(8);

  if(j == 0) {
    int size = state(1);
    // int size_diff = size - (1 + two_adults +  kids) ;
    double log_size_ratio = std::log(size) - std::log(1 + two_adults +  kids) ;

    double disp_inc = state(6) - state(2) - 2.5*sizes_sqm(size - 1); // income - rent
    disp_inc = disp_inc / sqrt(1 + two_adults + 0.5*kids) ; // equivalized (dep. on hh size)
    double hqual = state(0) * std::exp(-0.038 * state(3));

    // Utility function values
    mapped(0) += (disp_inc < .1) ?  -2.302585 * sw : std::log(disp_inc) * sw ; // income

    mapped(1) += hqual * sw; // hqual
    mapped(2) += sw * log_size_ratio ;
    mapped(3) += state(3) * sw ;// stay length

  } else {
    // moving costs

    mapped(nvars_u_c) += sw; // move dummy
    mapped(nvars_u_c + 1) += state(5) * sw; // variable moving costs (age)
    mapped(nvars_u_c + 2) += std::pow(state(5),2) * sw; // variable moving costs (age^2)
    mapped(nvars_u_c + 3) += two_adults * sw; // two adults
    mapped(nvars_u_c + 4) += kids * sw; // number of kids

    if(j <= 40) { // resets stay length to zero + rent & size & htype adjust

      int size = j0 / 10 + 1; // size
      // int size_diff = size - (1 + two_adults +  kids) ;
      double log_size_ratio = std::log(size) - std::log(1 + two_adults +  kids) ;

      int hqual = j - (size - 1)*10; // htype
      // double rent = sizes_sqm(size - 1) * state(8 + hqual);
      double rent = sizes_sqm(size - 1) * rent_hs(size - 1, hqual - 1); // rent

      double disp_inc = state(6) - rent - 2.5*sizes_sqm(size - 1); // income - rent
      disp_inc = disp_inc / sqrt(1 + two_adults + 0.5*kids) ; // equivalized (dep. on hh size)

      // Utility function values
      mapped(0) += (disp_inc < .1) ?  -2.302585 * sw : std::log(disp_inc) * sw ; // income

      mapped(1) += hqual * sw; // hqual

      mapped(2) += sw * log_size_ratio ;
      // mapped(3) += state(3) * sw // stay length = 0

    } else {

      int select = nvars_u_c + nvars_u_s + (j - 41) * nvars_v ;
      mapped(select + 0) += sw; // constant
      mapped(select + 1) += std::log(state(6)) * sw; // log income
      mapped(select + 2) += state(4)*1e-5 * sw; // wealth
      mapped(select + 3) += std::pow(state(4)*1e-5,2) * sw; // wealth^2
      mapped(select + 4) += state(5) * sw ; // age
      mapped(select + 5) += std::pow(state(5),2) * sw ; // age^2
      mapped(select + 6) += two_adults * sw; // two adults
      mapped(select + 7) += kids * sw; // number of kids

    }
  }

  return ;
}



// [[Rcpp::export]]
Rcpp::NumericVector f_x( // transition probabilities combined
    const int& M,
    const int& K,
    const int& E,
    Rcpp::NumericVector& state,
    Rcpp::NumericMatrix& adu_transition_coefs,
    Rcpp::NumericMatrix& kid_transition_coefs
) {
  // transition probabilities (adults, kids)
  Rcpp::NumericVector w_adu(M);
  Rcpp::NumericVector w_kid(K);
  // This does not depend on the choice.
  adu_transition_prob(w_adu, adu_transition_coefs, state); // prob: adults in t+1
  kid_transition_prob(w_kid, kid_transition_coefs, state);// prob: kids in t+1
  Rcpp::NumericVector f_x (E * M * K); // prob weights (do not depend on the choice)
  for(int m = 0; m < M; m++) {
    for(int k = 0; k < K; k++) {
      for(int e = 0; e < E; e++) {
        f_x(m * K * E + k * E + e) = w_adu(m) * w_kid(k) / E;
      }
    }
  }
  return f_x ;
}

//  path of different possible states
void get_paths(
    Rcpp::NumericVector& mapped_u1,
    Rcpp::NumericMatrix& States_j1, // (store all states -> psi)
    Rcpp::NumericMatrix& States_01, // (store all states -> psi)
    Rcpp::NumericVector& w_adu,
    Rcpp::NumericVector& w_kid,
    Rcpp::NumericVector& state,
    const int& j,
    const int& terminal_j,
    const int& E,
    const int& M,
    const int& K,
    const int& S,
    Rcpp::NumericVector& inc_transition_coefs,
    Rcpp::NumericVector& inc_transition_errors,
    Rcpp::NumericVector& sav_transition_coefs,
    Rcpp::NumericVector& sav_transition_errors,
    const int& nvars,
    const int& nvars_u_c,
    const int& nvars_u_s,
    const int& nvars_v,
    Rcpp::NumericVector& sizes_sqm,
    Rcpp::NumericMatrix& rent_hs_new,
    const bool& skip_hs,
    const double& g
) {
  // choice in period t has direct effects on period t util
  //   (moves are realized in period t already)

  int nstate = States_j1.ncol(); // number of state variables

  Rcpp::NumericVector state_01 (nstate); // deterministic state in t+1 if choice 0
  if(skip_hs) {
    update_det_state(state_01, state, -1, sizes_sqm, g);
  } else {
    update_det_state_hs(state_01, state, -1, sizes_sqm, rent_hs_new, g);
  }


  if(j < 40) {
    Rcpp::NumericVector state_j1 (nstate); // deterministic state in t+1 if choice j
    if(skip_hs) {
      update_det_state(state_j1, state, j, sizes_sqm, g);
    } else {
      update_det_state_hs(state_j1, state, j, sizes_sqm, rent_hs_new, g);
    }


    // adults m
    for(int m = 0; m < M; m++) {

      // income transition (depends on adults lead)
      Rcpp::NumericVector inc_lead (E);
      inc_transition(
        inc_lead,
        inc_transition_coefs, inc_transition_errors,
        state, m
      );

      // savings transition (depends on income lead & choice)
      Rcpp::NumericVector sav_lead_0 (E);
      sav_transition(
        sav_lead_0,
        sav_transition_coefs, sav_transition_errors,
        state, inc_lead,
        -1,
        sizes_sqm
      );
      Rcpp::NumericVector sav_lead_j (E);
      sav_transition(
        sav_lead_j,
        sav_transition_coefs, sav_transition_errors,
        state_j1, inc_lead, // previously: state, inc_lead
        j,
        sizes_sqm
      );

      state_01(7) = m; // two adults
      state_j1(7) = m; // two adults

      // kids k
      for(int k = 0; k < K; k++) {

        state_01(8) = k; // kids
        state_j1(8) = k; // kids


        double sw = w_adu(m) * w_kid(k) / E;

        for(int e = 0; e < E; e++) {
          state_01(6) = inc_lead(e); // income
          state_j1(6) = inc_lead(e); // income

          state_01(4) = sav_lead_0(e); // wealth / sav_lead already depends on state(4)
          state_j1(4) = sav_lead_j(e); // wealth / sav_lead already depends on state(4)

          // Utiility t+1
          //   baseline choice
          util(
            mapped_u1, state_01,
            terminal_j, S, -sw,
            nvars_u_c, nvars_u_s, nvars_v,
            sizes_sqm
          );
          //   alternative choice
          util(
            mapped_u1, state_j1,
            terminal_j, S, sw,
            nvars_u_c, nvars_u_s, nvars_v,
            sizes_sqm
          );

          // Store states
          int index_0 = m * K * E + k * E + e;
          int index_j = j * M * K * E + index_0 ;
          for(int i = 0; i < nstate; i++) {
            States_j1(index_j,i) = state_j1(i);
          }
          if(j == 0) {
            for(int i = 0; i < nstate; i++) {
              States_01(index_0,i) = state_01(i);
            }
          }
        }

      }
    }
  } else {

    for(int m = 0; m < M; m++) {

      // income transition (depends on adults lead)
      Rcpp::NumericVector inc_lead (E);
      inc_transition(
        inc_lead,
        inc_transition_coefs, inc_transition_errors,
        state, m
      );

      // savings transition (depends on income lead & choice)
      Rcpp::NumericVector sav_lead_0 (E);
      sav_transition(
        sav_lead_0,
        sav_transition_coefs, sav_transition_errors,
        state, inc_lead,
        -1,
        sizes_sqm
      );

      state_01(7) = m; // two adults

      // kids k
      for(int k = 0; k < K; k++) {

        state_01(8) = k; // kids

        double sw = w_adu(m) * w_kid(k) / E;

        for(int e = 0; e < E; e++) {
          state_01(6) = inc_lead(e); // income

          state_01(4) = sav_lead_0(e); // wealth / sav_lead already depends on state(4)

          // Utiility t+1
          //   baseline choice
          util(
            mapped_u1, state_01,
            terminal_j, S, -sw,
            nvars_u_c, nvars_u_s, nvars_v,
            sizes_sqm
          );

        }

      }
    }

  }
  //
  return ;
}


// store overall expected value of choice j
void mapped_vj(
    Rcpp::NumericMatrix& v,
    Rcpp::NumericMatrix& States_j1,
    Rcpp::NumericMatrix& States_01,
    Rcpp::NumericVector& f_x,
    Rcpp::NumericVector& w_adu,
    Rcpp::NumericVector& w_kid,
    Rcpp::NumericVector& state,
    const double& beta,
    const int& j,
    const int& terminal_j,
    const int& J,
    const int& E,
    const int& M,
    const int& K,
    const int& S,
    Rcpp::NumericVector& inc_transition_coefs,
    Rcpp::NumericVector& inc_transition_errors,
    Rcpp::NumericVector& sav_transition_coefs,
    Rcpp::NumericVector& sav_transition_errors,
    const int& nvars,
    const int& nvars_u_c,
    const int& nvars_u_s,
    const int& nvars_v,
    Rcpp::NumericVector& sizes_sqm,
    Rcpp::NumericMatrix& rent_hs,
    Rcpp::NumericMatrix& rent_hs_new,
    const bool& skip_hs,
    const double& g
) {
  Rcpp::NumericVector mapped_u0 (nvars); // t
  Rcpp::NumericVector mapped_u1 (nvars); // t+1


  // utility vector t mapped_u0
  //   baseline choice
  if(skip_hs) {
    util(
      mapped_u0, state,
      -1, S, -1,
      nvars_u_c, nvars_u_s, nvars_v,
      sizes_sqm
    );
    //   alternative choice
    util(
      mapped_u0, state,
      j, S, 1,
      nvars_u_c, nvars_u_s, nvars_v,
      sizes_sqm
    );
  } else {
    util_hs(
      mapped_u0, state,
      -1, S, -1,
      nvars_u_c, nvars_u_s, nvars_v,
      sizes_sqm,
      rent_hs
    );
    //   alternative choice
    util_hs(
      mapped_u0, state,
      j, S, 1,
      nvars_u_c, nvars_u_s, nvars_v,
      sizes_sqm,
      rent_hs
    );
  }


  // utility vector t+1 mapped_u1 + associated state
  get_paths(
    mapped_u1,
    States_j1, // attainable states in t+1 under alternative
    States_01, // attainable states in t+1 under baseline
    w_adu,
    w_kid,
    state,
    j,
    terminal_j,
    E,
    M,
    K,
    S,
    inc_transition_coefs,
    inc_transition_errors,
    sav_transition_coefs,
    sav_transition_errors,
    nvars,
    nvars_u_c,
    nvars_u_s,
    nvars_v,
    sizes_sqm,
    rent_hs_new,
    skip_hs,
    g
  );



  for(int s = 0; s < S; s++) {
    for(int i = 0; i < nvars_u_c; i++) {
      v(s*J + j, i) = mapped_u0(i) + beta * mapped_u1(i);
    }
    for(int i = 0; i < nvars_u_s; i++) {
      v(s*J + j, nvars_u_c + i) = mapped_u0(nvars_u_c + i) +
        beta * mapped_u1(nvars_u_c + i);
    }
    for(int i = 0; i < 3*nvars_v; i++) {
      v(s*J + j, nvars_u_c + nvars_u_s + i) = mapped_u0(nvars_u_c + nvars_u_s + i) +
        beta * mapped_u1(nvars_u_c + nvars_u_s + i);
    }
  }

  // nvars_u_c + nvars_u_s + (j - 41) * nvars_v
  // unobserved heterogeneity
  if(S == 8) {
    int select_s = nvars_u_c + nvars_u_s + 3*nvars_v - 1;
    int select_mc = nvars_u_c;
    int select_ex = nvars_u_c + nvars_u_s;
    int select_new = nvars_u_c + nvars_u_s + nvars_v;

    // s == 2 // moving cost dummy
    v(1*J + j, select_s + 1) = mapped_u0(select_mc) + beta * mapped_u1(select_mc); // moving cost dummy

    // s == 3 // existing dummy
    v(2*J + j, select_s + 2) = mapped_u0(select_ex) + beta * mapped_u1(select_ex); // existing dummy

    // s == 4 // new dummy
    v(3*J + j, select_s + 3) = mapped_u0(select_new) + beta * mapped_u1(select_new); // new dummy

    // s == 5 // moving cost dummy + existing dummy
    v(4*J + j, select_s + 1) = mapped_u0(select_mc) + beta * mapped_u1(select_mc); // moving cost dummy
    v(4*J + j, select_s + 2) = mapped_u0(select_ex) + beta * mapped_u1(select_ex); // existing dummy

    // s == 6 // moving cost dummy + new dummy
    v(5*J + j, select_s + 1) = mapped_u0(select_mc) + beta * mapped_u1(select_mc); // moving cost dummy
    v(5*J + j, select_s + 3) = mapped_u0(select_new) + beta * mapped_u1(select_new); // new dummy

    // s == 7 // existing dummy + new dummy
    v(6*J + j, select_s + 2) = mapped_u0(select_ex) + beta * mapped_u1(select_ex); // moving cost dummy
    v(6*J + j, select_s + 3) = mapped_u0(select_new) + beta * mapped_u1(select_new); // new dummy

    // s == 8 // moving cost dummy + existing dummy + new dummy
    v(7*J + j, select_s + 1) = mapped_u0(select_mc) + beta * mapped_u1(select_mc); // moving cost dummy
    v(7*J + j, select_s + 2) = mapped_u0(select_ex) + beta * mapped_u1(select_ex); // existing dummy
    v(7*J + j, select_s + 3) = mapped_u0(select_new) + beta * mapped_u1(select_new); // new dummy
  }
  return ;
}


// A & M (2011) Control term - choice path not optimal
void get_psi(
    Rcpp::NumericMatrix& psi,
    Rcpp::NumericMatrix& States_j1,
    Rcpp::NumericMatrix& States_01,
    Rcpp::NumericVector& f_x, // transition probabilities
    const int& S,
    const int& E,
    const int& K,
    const int& J,
    const int& M,
    const int& hid,
    const int& year,
    const bool& save_W
) {
  // compute log ccps
  Rcpp::NumericMatrix log_ccps (States_j1.nrow() + States_01.nrow(), S);
  Rcpp::Function f("compute_log_ccps"); // call to R function (uses updated q from workspace)
  log_ccps = f(States_j1, States_01, hid, year, f_x, save_W); // fill choice j first, baseline choice second

  int index_j, index_0 ;
  for(int j = 0; j < (J-3); j++) {
    for(int m = 0; m < M; m++) { // adults
      for(int k = 0; k < K; k++) { // kids
        for(int e = 0; e < E; e++) { // income x savings
          index_0 = m * K * E + k * E + e ;
          index_j = j * M * K * E + index_0;
          for(int s = 0; s < S; s++) {
            psi(j,s) -=  log_ccps(index_j,s) * f_x(index_0); // alternative j psi
          }
        }
      }
    }
  }
  index_j = (J - 3) * M * K * E ;
  for(int m = 0; m < M; m++) { // adults
    for(int k = 0; k < K; k++) { // kids
      for(int e = 0; e < E; e++) { // income x savings
        index_0 = m * K * E + k * E + e ;
        for(int s = 0; s < S; s++) {
          double tmp_psi = log_ccps(index_j + index_0,s) * f_x(index_0);
          for(int j = 0; j < J; j++) {
            psi(j,s) += tmp_psi ; // baseline choice psi
          }
        }
      }
    }
  }

  return;
}

// Same, but for simulation
void get_psi_sim(
    Rcpp::NumericMatrix& psi,
    Rcpp::NumericMatrix& States_j1,
    Rcpp::NumericMatrix& States_01,
    Rcpp::NumericVector& f_x, // transition probabilities
    const int& S,
    const int& E,
    const int& K,
    const int& J,
    const int& M
) {
  // compute log ccps
  Rcpp::NumericMatrix log_ccps (States_j1.nrow() + States_01.nrow(), 1);

  Rcpp::Function f("compute_log_ccps_sim"); // call to R function (uses updated q from workspace)
  log_ccps = f(States_j1, States_01); // fill choice j first, baseline choice second

  int index_j, index_0 ;
  for(int j = 0; j < (J-3); j++) {
    for(int m = 0; m < M; m++) { // adults
      for(int k = 0; k < K; k++) { // kids
        for(int e = 0; e < E; e++) { // income x savings
          index_0 = m * K * E + k * E + e ;
          index_j = j * M * K * E + index_0;
          for(int s = 0; s < S; s++) {
            psi(j,s) -=  log_ccps(index_j,0) * f_x(index_0); // alternative j psi
          }
        }
      }
    }
  }
  index_j = (J - 3) * M * K * E ;
  for(int m = 0; m < M; m++) { // adults
    for(int k = 0; k < K; k++) { // kids
      for(int e = 0; e < E; e++) { // income x savings
        index_0 = m * K * E + k * E + e ;
        for(int s = 0; s < S; s++) {
          double tmp_psi = log_ccps(index_j + index_0,0) * f_x(index_0);
          for(int j = 0; j < J; j++) {
            psi(j,s) += tmp_psi ; // baseline choice psi
          }
        }
      }
    }
  }

  return;
}


// Wrapper function: Compute valuations for each choice

// [[Rcpp::export]]
Rcpp::NumericVector mapped_v(
    Rcpp::NumericVector& state,
    const double& beta,
    const int& J,
    const int& terminal_j0,
    const int& nvars_u_c,
    const int& nvars_u_s,
    const int& nvars_v,
    const bool& compute_psi,
    Rcpp::NumericVector& inc_transition_coefs,
    Rcpp::NumericVector& inc_transition_errors,
    Rcpp::NumericVector& sav_transition_coefs,
    Rcpp::NumericVector& sav_transition_errors,
    Rcpp::NumericMatrix& adu_transition_coefs,
    Rcpp::NumericMatrix& kid_transition_coefs,
    const int& hid,
    const int& year,
    const bool& save_W,
    const bool& simulation,
    const int& s0,
    const int& S,
    Rcpp::NumericVector& sizes_sqm,
    Rcpp::NumericMatrix& rent_hs,
    Rcpp::NumericMatrix& rent_hs_new,
    const double& g
) {
  const bool skip_hs = (rent_hs.ncol() == 1) ;

  int nstate = state.length();

  const int E = sav_transition_errors.length(); // error pairs from income and savings processes
  const int M = 2; // # of states adults variable
  const int K = 3; // # of states kids variable

  const int nvars = nvars_u_c + nvars_u_s + nvars_v * 3;

  const int terminal_j = terminal_j0 - 1;

  const int S0 = (S == 8) ? 3 : 0 ;
  Rcpp::NumericMatrix v (J * S, nvars + 2 + S0);

  Rcpp::NumericMatrix States_j1 ((J-3) * E * M * K, nstate);
  Rcpp::NumericMatrix States_01 (E * M * K, nstate);

  // transition probabilities (adults, kids)
  Rcpp::NumericVector w_adu(M);
  Rcpp::NumericVector w_kid(K);
  // This does not depend on the choice.
  adu_transition_prob(w_adu, adu_transition_coefs, state); // prob: adults in t+1
  kid_transition_prob(w_kid, kid_transition_coefs, state);// prob: kids in t+1
  Rcpp::NumericVector f_x (E * M * K); // prob weights (do not depend on the choice)
  for(int m = 0; m < M; m++) {
    for(int k = 0; k < K; k++) {
      for(int e = 0; e < E; e++) {
        f_x(m * K * E + k * E + e) = w_adu(m) * w_kid(k) / E;
      }
    }
  }

  for(int j = 0; j < J; j++) {
    // initialize containers
    // Rcpp::NumericVector state_j1 (nstate); // deterministic state in t+1 if choice j
    // Rcpp::NumericVector state_01 (nstate); // deterministic state in t+1 if choice 0

    mapped_vj(
      v, // utility function vars (container)
      States_j1, // states in t+1 if alternative choice in t (container -> psi)
      States_01, // states in t+1 if baseline choice in t (container -> psi)
      f_x, // container for prob weights
      w_adu, // weights adults
      w_kid, // weights kids
      state, // current state
      beta, // discount factor
      j, // alternative j to be compared against baseline choice (j0 = -1)
      terminal_j, // terminal choice in t+1
      J, // total number of choices - 1
      E, // number of potential income x savings states
      M, // number of potential adult states
      K, // number of potential kids states
      S, // number of unobserved states
      inc_transition_coefs, // coefficients for income transition function
      inc_transition_errors, // errors income transition function
      sav_transition_coefs, // coefficients for savings transition function
      sav_transition_errors, // errors savings transition function
      nvars,
      nvars_u_c,
      nvars_u_s,
      nvars_v,
      sizes_sqm,
      rent_hs,
      rent_hs_new,
      skip_hs,
      g
    );

  }

  if(compute_psi) {
    Rcpp::NumericMatrix psi(J,S);
    if(simulation) {
      get_psi_sim(
        psi,
        States_j1,
        States_01,
        f_x,
        S,
        E,
        K,
        J,
        M
      );
    } else {
      get_psi(
        psi,
        States_j1,
        States_01,
        f_x,
        S,
        E,
        K,
        J,
        M,
        hid,
        year,
        save_W
      );
    }

    for(int s = 0; s < S; s++) {
      for(int j = 0; j < J; j++) {
        v(s*J + j, nvars + S0) = beta * psi(j,s);
      }
    }
  }


  for(int s = 0; s < S; s++) {
    for(int j = 0; j < J; j++) {
      v(s*J + j, nvars + 1 + S0) = s + 1;
    }
  }

  return v;
}




void get_paths_total(
    Rcpp::NumericVector& mapped_u1,
    Rcpp::NumericMatrix& States_j1, // (store all states -> psi)
    Rcpp::NumericMatrix& States_01, // (store all states -> psi)
    Rcpp::NumericVector& w_adu,
    Rcpp::NumericVector& w_kid,
    Rcpp::NumericVector& state,
    const int& j,
    const int& terminal_j,
    const int& E,
    const int& M,
    const int& K,
    const int& S,
    Rcpp::NumericVector& inc_transition_coefs,
    Rcpp::NumericVector& inc_transition_errors,
    Rcpp::NumericVector& sav_transition_coefs,
    Rcpp::NumericVector& sav_transition_errors,
    const int& nvars,
    const int& nvars_u_c,
    const int& nvars_u_s,
    const int& nvars_v,
    Rcpp::NumericVector& sizes_sqm,
    Rcpp::NumericMatrix& rent_hs_new,
    const bool& skip_hs,
    const double& g
) {
  // choice in period t has direct effects on period t util
  //   (moves are realized in period t already)

  int nstate = States_j1.ncol(); // number of state variables

  Rcpp::NumericVector state_01 (nstate); // deterministic state in t+1 if choice 0
  if(skip_hs) {
    update_det_state(state_01, state, -1, sizes_sqm, g);
  } else {
    update_det_state_hs(state_01, state, -1, sizes_sqm, rent_hs_new, g);
  }


  if(j < 40) {
    Rcpp::NumericVector state_j1 (nstate); // deterministic state in t+1 if choice j
    if(skip_hs) {
      update_det_state(state_j1, state, j, sizes_sqm, g);
    } else {
      update_det_state_hs(state_j1, state, j, sizes_sqm, rent_hs_new, g);
    }


    // adults m
    for(int m = 0; m < M; m++) {

      // income transition (depends on adults lead)
      Rcpp::NumericVector inc_lead (E);
      inc_transition(
        inc_lead,
        inc_transition_coefs, inc_transition_errors,
        state, m
      );

      // savings transition (depends on income lead & choice)
      Rcpp::NumericVector sav_lead_0 (E);
      sav_transition(
        sav_lead_0,
        sav_transition_coefs, sav_transition_errors,
        state, inc_lead,
        -1,
        sizes_sqm
      );
      Rcpp::NumericVector sav_lead_j (E);
      sav_transition(
        sav_lead_j,
        sav_transition_coefs, sav_transition_errors,
        state_j1, inc_lead, // previously: state, inc_lead
        j,
        sizes_sqm
      );

      state_01(7) = m; // two adults
      state_j1(7) = m; // two adults

      // kids k
      for(int k = 0; k < K; k++) {

        state_01(8) = k; // kids
        state_j1(8) = k; // kids

        // state_01(19) = ((k > 0) & (state(8) > 0)) ? state(19) + 1 : 0 ; // age oldest child
        // state_j1(19) = state_01(19); // age oldest child

        double sw = w_adu(m) * w_kid(k) / E;

        for(int e = 0; e < E; e++) {
          state_01(6) = inc_lead(e); // income
          state_j1(6) = inc_lead(e); // income

          state_01(4) = state(4) + sav_lead_0(e); // wealth
          state_j1(4) = state(4) + sav_lead_j(e); // wealth

          // Utiility t+1
          //   alternative choice
          util(
            mapped_u1, state_j1,
            terminal_j, S, sw,
            nvars_u_c, nvars_u_s, nvars_v,
            sizes_sqm
          );

          // Store states
          int index_0 = m * K * E + k * E + e;
          int index_j = j * M * K * E + index_0 ;
          for(int i = 0; i < nstate; i++) {
            States_j1(index_j,i) = state_j1(i);
          }
          if(j == 0) {
            for(int i = 0; i < nstate; i++) {
              States_01(index_0,i) = state_01(i);
            }
          }
        }

      }
    }
  }
  //
  return ;
}




void mapped_vj_total(
    Rcpp::NumericMatrix& v,
    Rcpp::NumericMatrix& States_j1,
    Rcpp::NumericMatrix& States_01,
    Rcpp::NumericVector& f_x,
    Rcpp::NumericVector& w_adu,
    Rcpp::NumericVector& w_kid,
    Rcpp::NumericVector& state,
    const double& beta,
    const int& j,
    const int& terminal_j,
    const int& J,
    const int& E,
    const int& M,
    const int& K,
    const int& S,
    Rcpp::NumericVector& inc_transition_coefs,
    Rcpp::NumericVector& inc_transition_errors,
    Rcpp::NumericVector& sav_transition_coefs,
    Rcpp::NumericVector& sav_transition_errors,
    const int& nvars,
    const int& nvars_u_c,
    const int& nvars_u_s,
    const int& nvars_v,
    Rcpp::NumericVector& sizes_sqm,
    Rcpp::NumericMatrix& rent_hs,
    Rcpp::NumericMatrix& rent_hs_new,
    const bool& skip_hs,
    const double& g
) {
  Rcpp::NumericVector mapped_u0 (nvars); // t
  Rcpp::NumericVector mapped_u1 (nvars); // t+1


  // utility vector t mapped_u0
  //   baseline choice
  if(skip_hs) {
    //   alternative choice
    util(
      mapped_u0, state,
      j, S, 1,
      nvars_u_c, nvars_u_s, nvars_v,
      sizes_sqm
    );
  } else {
    //   alternative choice
    util_hs(
      mapped_u0, state,
      j, S, 1,
      nvars_u_c, nvars_u_s, nvars_v,
      sizes_sqm,
      rent_hs
    );
  }


  // utility vector t+1 mapped_u1 + associated state
  get_paths_total(
    mapped_u1,
    States_j1, // attainable states in t+1 under alternative
    States_01, // attainable states in t+1 under baseline
    w_adu,
    w_kid,
    state,
    j,
    terminal_j,
    E,
    M,
    K,
    S,
    inc_transition_coefs,
    inc_transition_errors,
    sav_transition_coefs,
    sav_transition_errors,
    nvars,
    nvars_u_c,
    nvars_u_s,
    nvars_v,
    sizes_sqm,
    rent_hs_new,
    skip_hs,
    g
  );



  for(int s = 0; s < S; s++) {
    for(int i = 0; i < nvars_u_c; i++) {
      v(s*J + j, i) = mapped_u0(i) + beta * mapped_u1(i);
    }
    for(int i = 0; i < nvars_u_s; i++) {
      v(s*J + j, nvars_u_c + i) = mapped_u0(nvars_u_c + i) +
        beta * mapped_u1(nvars_u_c + i);
    }
    for(int i = 0; i < 3*nvars_v; i++) {
      v(s*J + j, nvars_u_c + nvars_u_s + i) = mapped_u0(nvars_u_c + nvars_u_s + i) +
        beta * mapped_u1(nvars_u_c + nvars_u_s + i);
    }
  }

  // nvars_u_c + nvars_u_s + (j - 41) * nvars_v
  // unobserved heterogeneity
  if(S == 8) {
    int select_s = nvars_u_c + nvars_u_s + 3*nvars_v - 1;
    int select_mc = nvars_u_c;
    int select_ex = nvars_u_c + nvars_u_s;
    int select_new = nvars_u_c + nvars_u_s + nvars_v;

    // s == 2 // moving cost dummy
    v(1*J + j, select_s + 1) = mapped_u0(select_mc) + beta * mapped_u1(select_mc); // moving cost dummy

    // s == 3 // existing dummy
    v(2*J + j, select_s + 2) = mapped_u0(select_ex) + beta * mapped_u1(select_ex); // existing dummy

    // s == 4 // new dummy
    v(3*J + j, select_s + 3) = mapped_u0(select_new) + beta * mapped_u1(select_new); // new dummy

    // s == 5 // moving cost dummy + existing dummy
    v(4*J + j, select_s + 1) = mapped_u0(select_mc) + beta * mapped_u1(select_mc); // moving cost dummy
    v(4*J + j, select_s + 2) = mapped_u0(select_ex) + beta * mapped_u1(select_ex); // existing dummy

    // s == 6 // moving cost dummy + new dummy
    v(5*J + j, select_s + 1) = mapped_u0(select_mc) + beta * mapped_u1(select_mc); // moving cost dummy
    v(5*J + j, select_s + 3) = mapped_u0(select_new) + beta * mapped_u1(select_new); // new dummy

    // s == 7 // existing dummy + new dummy
    v(6*J + j, select_s + 2) = mapped_u0(select_ex) + beta * mapped_u1(select_ex); // moving cost dummy
    v(6*J + j, select_s + 3) = mapped_u0(select_new) + beta * mapped_u1(select_new); // new dummy

    // s == 8 // moving cost dummy + existing dummy + new dummy
    v(7*J + j, select_s + 1) = mapped_u0(select_mc) + beta * mapped_u1(select_mc); // moving cost dummy
    v(7*J + j, select_s + 2) = mapped_u0(select_ex) + beta * mapped_u1(select_ex); // existing dummy
    v(7*J + j, select_s + 3) = mapped_u0(select_new) + beta * mapped_u1(select_new); // new dummy
  }
  return ;
}



// [[Rcpp::export]]
Rcpp::NumericVector mapped_v_total(
    Rcpp::NumericVector& state,
    const double& beta,
    const int& J,
    const int& terminal_j0,
    const int& nvars_u_c,
    const int& nvars_u_s,
    const int& nvars_v,
    const bool& compute_psi,
    Rcpp::NumericVector& inc_transition_coefs,
    Rcpp::NumericVector& inc_transition_errors,
    Rcpp::NumericVector& sav_transition_coefs,
    Rcpp::NumericVector& sav_transition_errors,
    Rcpp::NumericMatrix& adu_transition_coefs,
    Rcpp::NumericMatrix& kid_transition_coefs,
    const int& hid,
    const int& year,
    const bool& save_W,
    const bool& simulation,
    const int& s0,
    const int& S,
    Rcpp::NumericVector& sizes_sqm,
    Rcpp::NumericMatrix& rent_hs,
    Rcpp::NumericMatrix& rent_hs_new,
    const double& g
) {
  const bool skip_hs = (rent_hs.ncol() == 1) ;

  int nstate = state.length();

  const int E = sav_transition_errors.length(); // error pairs from income and savings processes
  const int M = 2; // # of states adults variable
  const int K = 3; // # of states kids variable

  const int nvars = nvars_u_c + nvars_u_s + nvars_v * 3;

  const int terminal_j = terminal_j0 - 1;

  const int S0 = (S == 8) ? 3 : 0 ;
  Rcpp::NumericMatrix v (J * S, nvars + 2 + S0);

  Rcpp::NumericMatrix States_j1 ((J-3) * E * M * K, nstate);
  Rcpp::NumericMatrix States_01 (E * M * K, nstate);

  // transition probabilities (adults, kids)
  Rcpp::NumericVector w_adu(M);
  Rcpp::NumericVector w_kid(K);
  // This does not depend on the choice.
  adu_transition_prob(w_adu, adu_transition_coefs, state); // prob: adults in t+1
  kid_transition_prob(w_kid, kid_transition_coefs, state);// prob: kids in t+1
  Rcpp::NumericVector f_x (E * M * K); // prob weights (do not depend on the choice)
  for(int m = 0; m < M; m++) {
    for(int k = 0; k < K; k++) {
      for(int e = 0; e < E; e++) {
        f_x(m * K * E + k * E + e) = w_adu(m) * w_kid(k) / E;
      }
    }
  }

  for(int j = 0; j < J; j++) {
    // initialize containers
    // Rcpp::NumericVector state_j1 (nstate); // deterministic state in t+1 if choice j
    // Rcpp::NumericVector state_01 (nstate); // deterministic state in t+1 if choice 0

    mapped_vj_total(
      v, // utility function vars (container)
      States_j1, // states in t+1 if alternative choice in t (container -> psi)
      States_01, // states in t+1 if baseline choice in t (container -> psi)
      f_x, // container for prob weights
      w_adu, // weights adults
      w_kid, // weights kids
      state, // current state
      beta, // discount factor
      j, // alternative j to be compared against baseline choice (j0 = -1)
      terminal_j, // terminal choice in t+1
      J, // total number of choices - 1
      E, // number of potential income x savings states
      M, // number of potential adult states
      K, // number of potential kids states
      S, // number of unobserved states
      inc_transition_coefs, // coefficients for income transition function
      inc_transition_errors, // errors income transition function
      sav_transition_coefs, // coefficients for savings transition function
      sav_transition_errors, // errors savings transition function
      nvars,
      nvars_u_c,
      nvars_u_s,
      nvars_v,
      sizes_sqm,
      rent_hs,
      rent_hs_new,
      skip_hs,
      g
    );

  }

  if(compute_psi) {
    Rcpp::NumericMatrix psi(J,S);
    if(simulation) {
      get_psi_sim(
        psi,
        States_j1,
        States_01,
        f_x,
        S,
        E,
        K,
        J,
        M
      );
    } else {
      get_psi(
        psi,
        States_j1,
        States_01,
        f_x,
        S,
        E,
        K,
        J,
        M,
        hid,
        year,
        save_W
      );
    }

    for(int s = 0; s < S; s++) {
      for(int j = 0; j < J; j++) {
        v(s*J + j, nvars + S0) = beta * psi(j,s);
      }
    }
  }


  for(int s = 0; s < S; s++) {
    for(int j = 0; j < J; j++) {
      v(s*J + j, nvars + 1 + S0) = s + 1;
    }
  }

  return v;
}







double Lik_Choice(
    Rcpp::NumericMatrix& x,
    Rcpp::NumericVector& theta,
    const int& j0
) {
  double V = 1;
  Rcpp::NumericVector v(x.nrow() + 1);
  v(0) = 1;
  for(int j = 1; j < v.length(); j++) {
    for(int i = 0; i < theta.length(); i++) {
      v(j) += theta(i) * x(j-1, i) ;
    }
    v(j) += x(j-1, theta.length()) ; // psi
    v(j) = std::exp(v(j));
    V += v(j);
  }
  double LL = v(j0)/V;
  return LL;
}



//  eq 5.3

// state vars
// > data.frame(pos = 0:(length(state.vars)-1), var = state.vars)
//      pos                  var
//   1    0                htype
//   2    1          stay_length
//   3    2          lock_in_inc
//   4    3               wealth
//   5    4                  age
//   6    5     income_class_ror
//   7    6 income_class_ror_lag
//   8    7           match_qual

// [[Rcpp::export]]
Rcpp::NumericVector update_qn(
    Rcpp::NumericVector& choices, // sequence of choices for HH n, t = 1,..,T
    Rcpp::NumericVector& exp_v,
    Rcpp::NumericVector& pi, // distribution of unobserved types in the data
    const int& J
) {
  int T = choices.length();
  const int S = pi.length();

  double L = 0;
  Rcpp::NumericVector Lns(S);
  for(int s = 0; s < S; s++) {
    double Ls = pi(s) ;
    Lns(s) = pi(s) ;
    for(int t = 0; t < T; t++) {
      int index = t*J*S + s*J;
      Rcpp::NumericVector v0 = exp_v[
        Rcpp::Range(index, index + J - 1)
      ];
      double V = Rcpp::sum(v0);
      V++;
      double Lt = (choices(t) > 0) ? v0(choices(t) - 1) / V : 1/V ;
      Ls *= Lt ;
      Lns(s) *= Lt;
    }
    L += Ls;
  }
  Rcpp::NumericVector qn(S);
  for(int s = 0; s < S; s++) {
    qn(s) = Lns(s) / L ;
  }
  return qn;
}
//
void Choice_Prob(
    Rcpp::NumericVector& prob,
    Rcpp::NumericMatrix& x,
    Rcpp::NumericVector& theta
) {
  double V = 1.0;
  prob(0) = 1;

  for(int j = 1; j < prob.length(); j++) {

    for(int i = 0; i < theta.length(); i++) {
      prob(j) += theta(i) * x(j - 1, i) ;
    }

    prob(j) += x(j - 1, theta.length()); // psi
    prob(j) = std::exp(prob(j));
    V += prob(j);
  }

  for(int j = 0; j < prob.length(); j++) {
    prob(j) = prob(j)/V;
  }

  return ;
}
//
void Gradient(
    Rcpp::NumericVector& grad,
    Rcpp::NumericVector& prob,
    const int& choice,
    Rcpp::NumericMatrix& v0
) {
  int J = v0.nrow();
  int K = grad.length();
  for(int k = 0; k < K; k++) {
    if(choice > 0) grad(k) += v0(choice - 1,k) ;
    for(int j = 0; j < J; j++) {
      grad(k) -= prob(j + 1) * v0(j,k);
    }
  }
  return ;
}
//

// [[Rcpp::export]]
double LogLik_contribution(
    Rcpp::NumericMatrix& v_n, // mapped utility vars for HH n's choices at t = 1,..,T, s=1,..,S
    Rcpp::NumericVector& theta,
    Rcpp::NumericMatrix& qn, // conditional distribution of unobserved types, hh n
    Rcpp::NumericVector& start_index, // hid start row numbers
    const int& J
) {
  // transitions don't matter since they
  // do not depend on s and can thus be estimated separately
  const int S = qn.ncol();
  const int H = start_index.length() - 1;

  double LL = 0.0;
  // Rcpp::NumericVector grad(theta.length()) ; // gradient
  // Rcpp::NumericMatrix H(theta.length(), theta.length()); // hessian

  for(int h = 0; h < H; h++) {
    const int T = (start_index[h+1] - start_index[h]) / (S * J);

    for(int s = 0; s < S; s++) {

      double L_s = 0.0 ;

      for(int t = 0; t < T; t++) {

        const int index = start_index[h] + t*J*S + s * J;
        Rcpp::NumericMatrix v0 = v_n(
          Rcpp::Range(index, index + J - 1),
          Rcpp::Range(4, v_n.ncol() - 1)
        );

        Rcpp::NumericVector prob(v0.nrow() + 1);
        Choice_Prob(prob, v0, theta);

        L_s += std::log( prob(v_n(index,2)) ); // v_n(index,2) is the actual choice

      }

      LL += qn(h,s) * L_s ;
    }
  }
  return LL;
}
//
// [[Rcpp::export]]
Rcpp::NumericVector gradient_contribution(
    Rcpp::NumericMatrix& v_n, // mapped utility vars for HH n's choices at t = 1,..,T, s=1,..,S
    Rcpp::NumericVector& theta,
    Rcpp::NumericMatrix& qn, // conditional distribution of unobserved types, hh n
    Rcpp::NumericVector& start_index, // hid start row numbers
    const int& J
) {
  // transitions don't matter since they
  // do not depend on s and can thus be estimated separately
  const int S = qn.ncol();
  const int H = start_index.length() - 1;

  Rcpp::NumericVector grad(theta.length()) ; // gradient


  for(int h = 0; h < H; h++) {
    const int T = (start_index[h+1] - start_index[h]) / (S * J);

    for(int s = 0; s < S; s++) {
      Rcpp::NumericVector grad_s(theta.length()) ; // gradient

      for(int t = 0; t < T; t++) {
        const int index = start_index[h] + t*J*S + s * J;
        Rcpp::NumericMatrix v0 = v_n(
          Rcpp::Range(index, index + J - 1),
          Rcpp::Range(4, v_n.ncol() - 1)
        );

        Rcpp::NumericVector prob(v0.nrow() + 1);
        Choice_Prob(prob, v0, theta);

        Gradient(grad_s, prob, v_n(index,2), v0);
      }

      for(int k = 0; k < grad.length(); k++) {
        grad(k) += qn(h, s) * grad_s(k);
      }
    }

  }

  return grad;
}
//
//

//
//

// Rent Gradient Demand / Supply equilibrium


// [[Rcpp::export]]
double p_42(
    Rcpp::NumericVector& state,
    Rcpp::NumericVector& coefs
) {
  Rcpp::NumericVector x (coefs.length());

  // These inputs depend on the estimated choice model
  //  used to parameterize p_42 in the simulations.
  //   choice 42 == choice to buy existing home
  // Coefficients:
  //   Estimate Std. Error z value Pr(>|z|)
  //   (Intercept)      -5.868e+00  4.126e-01 -14.221  < 2e-16 ***
  //   income            9.909e-04  1.629e-04   6.083 1.18e-09 ***
  //   htype_adj        -3.194e-02  2.442e-02  -1.308 0.190956
  //   size              3.552e-02  1.037e-01   0.343 0.731940
  //   stay_length       2.435e-01  7.273e-02   3.348 0.000815 ***
  //   age              -2.212e-02  5.966e-03  -3.707 0.000210 ***
  //   two_adults        4.673e-01  1.902e-01   2.456 0.014032 *
  //   kids              1.934e-01  8.121e-02   2.382 0.017215 *
  //   gwealth_imp       3.644e-06  1.077e-06   3.384 0.000714 ***
  //   rent             -1.136e-04  3.104e-04  -0.366 0.714269
  //   I(income^2)      -7.887e-08  1.818e-08  -4.338 1.44e-05 ***
  //   I(stay_length^2) -3.282e-02  9.846e-03  -3.333 0.000858 ***
  //   age:gwealth_imp  -5.021e-08  1.982e-08  -2.533 0.011310 *
  x(0) = 1 ; // intercept

  x(1) = state(6) ; // income

  x(2) = state(0) ; // htype_adj
  x(3) = state(1) ; // size
  x(4) = state(3) ; // stay_length

  x(5) = state(5) ; // age
  x(6) = state(7) ; // two adults
  x(7) = state(8) ; // kids

  x(8) = state(4) ;  // gwealth_imp

  x(9) = state(2) ; // rent

  x(10) = std::pow(state(6), 2) ; // I(income^2)
  x(11) = std::pow(state(3), 2) ; // I(stay_length^2)

  x(12) = state(5) * state(4) ; // age x gwealth_imp

  //
  double V = 0;
  for(int k = 0; k < x.length(); k++) {
    V += x(k) * coefs(k) ;
  }
  V = std::exp(V);
  double prob = V / (1 + V) ;
  prob = (prob > 0.0001) ? prob : 0.0001 ;
  return prob;
}


double g_df42_dx(
    Rcpp::NumericVector& state,
    Rcpp::NumericVector& coefs,
    const double& dwealth_dr,
    const int& h,
    Rcpp::NumericVector& sizes_sqm
) {
  // int size = h / 10; // size
  // int hqual = h - size*10 + 1; // htype

  double dx = ( // effect working through gwealth_imp
    coefs(8) + // main effect of wealth
      coefs(12) * state(5) // interaction with age
  ) * dwealth_dr ;

  dx +=  ( // effect working through rent0
    coefs(9)
  ) * sizes_sqm( state(1) - 1 ) ;

  return dx;
}

// derivative of current-period utility w.r.t. r
void du_dr(
    Rcpp::NumericMatrix& dv_dr,
    Rcpp::NumericVector& state,
    Rcpp::NumericVector& theta,
    const int& j0,
    Rcpp::NumericVector& sizes_sqm,
    Rcpp::NumericMatrix& rent_hs
) {
  int two_adults = state(7);
  int kids = state(8);
  // int sel = (kids == 0) ? 2 : 0;

  int size = j0 / 10 + 1; // size
  int hqual = j0 - (size - 1)*10 + 1; // htype
  double rent = sizes_sqm(size - 1) * rent_hs(size - 1, hqual - 1); // rent

  double disp_inc = state(6) - rent - 2.5*sizes_sqm(size - 1); // income - rent - costs
  double scale = 1 / std::sqrt(1 + two_adults + 0.5*kids); //
  // disp_inc = disp_inc * scale ;
  //
  // dv_dr(j0, j0) -= theta(sel) * sizes_sqm(size - 1) * scale ;
  // dv_dr(j0, j0) -= theta(sel+1) * (2 * disp_inc * scale) * sizes_sqm(size - 1) * scale ;

  // double inv_disp_inc = (disp_inc > .1) ? 1/disp_inc : 0 ;
  double d_du = (disp_inc * scale > .1) ?  -sizes_sqm(size - 1) / disp_inc : 0 ;
  dv_dr(j0, j0) += theta(0) * d_du ;

  return ;
}

void dv1_dr(
    Rcpp::NumericMatrix& dv_dr,
    Rcpp::NumericVector& state_0,
    Rcpp::NumericVector& state_j,
    Rcpp::NumericVector& theta,
    Rcpp::NumericVector& sav_transition_coefs,
    const double& sw,
    const int& j0,
    Rcpp::NumericVector& sizes_sqm,
    Rcpp::NumericVector& p42_coefs
) {
  if(j0 < 40) {
    // du42_dr (impact of rent change on future terminal utility)
    // Impact on utility due to wealth change in t+1 from rent change in t
    //  This term depends on wealth, but it does not depend on the rent level.
    //  The wealth transition depends on the rent level.
    double dwealth_dr = 12 * sizes_sqm( state_j(1) - 1 ) * sav_transition_coefs(2) ;
    double du42_dwealth = theta(11) * 1e-5 + theta(12) * state_j(4) * 2e-10 ;
    double du42_dr =  du42_dwealth * dwealth_dr;

    // Impact on correction term under alternative
    //  g_df42_dx() depends on wealth.
    //  g_df42_dx() may also depend on the rent, but currently doesn't.
    //    (Try to include the rent if it works w/o the rent.)
    double p42_alt = p_42(state_j, p42_coefs) ;// probability to choose 42
    double dp42_dr_alt = (
      p42_alt > 0.0001
    ) ? (1 - p42_alt) * g_df42_dx(state_j, p42_coefs, dwealth_dr, j0, sizes_sqm) : 0;
    // p = exp(f(x)) / (1 + exp(f(x)))
    // dlog(p)/dx = 1/p * dp/dx = (1/p) * (1 - p)^2 * exp(f(x)) * f'(x) = (1 - p) * f'(x)
    dv_dr(j0, j0) += sw * (du42_dr - dp42_dr_alt);
  }

  return;
}



void dvj_drk(
    Rcpp::NumericMatrix& dv_dr,
    Rcpp::NumericVector& theta,
    Rcpp::NumericVector& w_adu,
    Rcpp::NumericVector& w_kid,
    Rcpp::NumericVector& state,
    const double& beta,
    const int& j,
    const int& terminal_j,
    const int& E,
    const int& M,
    const int& K,
    const int& S,
    Rcpp::NumericVector& inc_transition_coefs,
    Rcpp::NumericVector& inc_transition_errors,
    Rcpp::NumericVector& sav_transition_coefs,
    Rcpp::NumericVector& sav_transition_errors,
    Rcpp::NumericVector& p42_coefs,
    const int& nvars,
    const int& nvars_u_c,
    const int& nvars_u_s,
    const int& nvars_v,
    Rcpp::NumericVector& sizes_sqm,
    Rcpp::NumericMatrix& rent_hs,
    Rcpp::NumericMatrix& rent_hs_new,
    const double& g
) {

  int nstate = state.length(); // number of state variables

  Rcpp::NumericVector state_01 (nstate); // deterministic state in t+1 if choice 0
  update_det_state_hs(state_01, state, -1, sizes_sqm, rent_hs_new, g);

  Rcpp::NumericVector state_j1 (nstate); // deterministic state in t+1 if choice j


  if(j < 40) { // for j >= 40, there still is an effect of dr on p_42 under the baseline choice
    du_dr(
      dv_dr,
      state, theta, j,
      sizes_sqm, rent_hs
    );

    update_det_state_hs(state_j1, state, j, sizes_sqm, rent_hs_new, g);

    // adults m
    for(int m = 0; m < M; m++) {

      // income transition (depends on adults lead)
      Rcpp::NumericVector inc_lead (E);
      inc_transition(
        inc_lead,
        inc_transition_coefs, inc_transition_errors,
        state, m
      );

      // savings transition (depends on income lead & choice)
      Rcpp::NumericVector sav_lead_0 (E);
      sav_transition(
        sav_lead_0,
        sav_transition_coefs, sav_transition_errors,
        state, inc_lead,
        -1,
        sizes_sqm
      );
      Rcpp::NumericVector sav_lead_j (E);
      sav_transition(
        sav_lead_j,
        sav_transition_coefs, sav_transition_errors,
        state_j1, inc_lead,
        j,
        sizes_sqm
      );

      state_01(7) = m; // two adults
      state_j1(7) = m; // two adults

      // kids k
      for(int k = 0; k < K; k++) {

        state_01(8) = k; // kids
        state_j1(8) = k; // kids

        // state_01(19) = ((k > 0) & (state(8) > 0)) ? state(19) + 1 : 0 ; // age oldest child
        // state_j1(19) = state_01(19); // age oldest child

        double sw = beta * w_adu(m) * w_kid(k) / E;

        for(int e = 0; e < E; e++) {
          state_01(6) = inc_lead(e); // income
          state_j1(6) = inc_lead(e); // income

          state_01(4) = sav_lead_0(e); // wealth
          state_j1(4) = sav_lead_j(e); // wealth

          // dv_(t+1)/dr
          dv1_dr(
              dv_dr,
              state_01, state_j1,
              theta, sav_transition_coefs,
              sw, j,
              sizes_sqm,
              p42_coefs
          );

        }

      }
    }
  }
  // else {
  //
  //   for(int m = 0; m < M; m++) {
  //
  //     // income transition (depends on adults lead)
  //     Rcpp::NumericVector inc_lead (E);
  //     inc_transition(
  //       inc_lead,
  //       inc_transition_coefs, inc_transition_errors,
  //       state, m
  //     );
  //
  //     // savings transition (depends on income lead & choice)
  //     Rcpp::NumericVector sav_lead_0 (E);
  //     sav_transition(
  //       sav_lead_0,
  //       sav_transition_coefs, sav_transition_errors,
  //       state, inc_lead,
  //       -1,
  //       sizes_sqm
  //     );
  //
  //     state_01(7) = m; // two adults
  //
  //     // kids k
  //     for(int k = 0; k < K; k++) {
  //
  //       state_01(8) = k; // kids
  //
  //       double sw = beta * w_adu(m) * w_kid(k) / E;
  //
  //       for(int e = 0; e < E; e++) {
  //         state_01(6) = inc_lead(e); // income
  //
  //         state_01(4) = sav_lead_0(e); // wealth
  //
  //         // dv_(t+1)/dr
  //         dv1_dr(
  //             dv_dr,
  //             state_01, state_j1,
  //             theta, sav_transition_coefs,
  //             sw, j,
  //             sizes_sqm,
  //             p42_coefs
  //         );
  //       }
  //
  //     }
  //   }
  //
  // }
  //
  return ;

}



// [[Rcpp::export]]
Rcpp::NumericMatrix dv_ccp(
    Rcpp::NumericVector& state,
    const double& beta,
    const int& J, // total number of choices - 1 (excl. baseline, stay choice)
    const int& terminal_j0,
    const int& nvars_u_c,
    const int& nvars_u_s,
    const int& nvars_v,
    const bool& compute_psi,
    Rcpp::NumericVector& inc_transition_coefs,
    Rcpp::NumericVector& inc_transition_errors,
    Rcpp::NumericVector& sav_transition_coefs,
    Rcpp::NumericVector& sav_transition_errors,
    Rcpp::NumericMatrix& adu_transition_coefs,
    Rcpp::NumericMatrix& kid_transition_coefs,
    Rcpp::NumericVector& p42_coefs,
    const int& hid,
    const int& year,
    const bool& save_W,
    const bool& simulation,
    const int& s0,
    const int& S,
    Rcpp::NumericVector& sizes_sqm,
    Rcpp::NumericMatrix& rent_hs,
    Rcpp::NumericMatrix& rent_hs_new,
    Rcpp::NumericVector& theta,
    Rcpp::NumericVector& pn,
    const double& g
) {


  const int E = sav_transition_errors.length(); // error pairs from income and savings processes
  const int M = 2; // # of states adults variable
  const int K = 3; // # of states kids variable

  const int nvars = nvars_u_c + nvars_u_s * S + nvars_v * 3;

  const int terminal_j = terminal_j0 - 1;

  // Container: Jacobi Matrix dv/dr
  Rcpp::NumericMatrix dv_dr (J, 40);

  // transition probabilities (adults, kids)
  Rcpp::NumericVector w_adu(M);
  Rcpp::NumericVector w_kid(K);
  // This does not depend on the choice.
  adu_transition_prob(w_adu, adu_transition_coefs, state); // prob: adults in t+1
  kid_transition_prob(w_kid, kid_transition_coefs, state);// prob: kids in t+1

  for(int j = 0; j < J; j++) {
    dvj_drk(
      dv_dr, // container for dv/dr
      theta,  // estimated parameter (theta)
      w_adu, // weights adults
      w_kid, // weights kids
      state, // current state
      beta, // discount factor
      j, // alternative j to be compared against baseline choice (j0 = -1)
      terminal_j, // terminal choice in t+1
      E, // number of potential income x savings states
      M, // number of potential adult states
      K, // number of potential kids states
      S, // number of unobserved states
      inc_transition_coefs, // coefficients for income transition function
      inc_transition_errors, // errors income transition function
      sav_transition_coefs, // coefficients for savings transition function
      sav_transition_errors, // errors savings transition function
      p42_coefs,
      nvars,
      nvars_u_c,
      nvars_u_s,
      nvars_v,
      sizes_sqm,
      rent_hs,
      rent_hs_new,
      g
    );
  }


  Rcpp::NumericMatrix d_pn (pn.length(), dv_dr.ncol());
  for(int h = 0; h < 40; h++) {
    double netsum = 0 ;
    for(int j = 0; j < (J * S); j++) {
      netsum += pn(j+1) * dv_dr(j,h);
    }
    for(int j = 0; j < (J * S); j++) {
      d_pn(j+1,h) = pn(j+1) * (dv_dr(j,h) - netsum);
    }
    d_pn(0,h) = -pn(0) * netsum;
  }

  //return dv_dr;
  return d_pn;
}












// /////////////// Variance-Covariance Matrix ///////////////
//
// // [[Rcpp::export]]
// Rcpp::List vcov_parts(
//   Rcpp::NumericVector& choices_n, // sequence of choices for HH n, t = 1,..,T
//   Rcpp::NumericMatrix& states_n, // sequence of states for HH n, t = 1,..,T+1
//   Rcpp::NumericMatrix& v_n, // mapped utility vars for HH n's choices at t = 1,..,T, s=1,..,S
//   Rcpp::NumericVector& pdiff_n, // discounted diff of inverse ccp in t+1, t = 1,..,T
//   Rcpp::NumericVector& w_res_n, // empirical ccp choice residuals * HH weight in CCP estimator
//   Rcpp::NumericVector& q_n, // distribution of unobserved types for HH n
//   double& w_n, // avg. weight of HH n in the CCP estimator
//   Rcpp::NumericVector& theta, // utility function parameters (w/o psi)
//   Rcpp::NumericVector& pi, // distribution of unobserved types in the data
//   Rcpp::NumericVector& inc_transition_coefs,
//   Rcpp::NumericVector& match_transition_coefs,
//   Rcpp::NumericVector& savings_transition_coefs,
//   const int& J
// ) {
//   const int T = choices_n.length();
//   const int E = 10;
//   const int M = 3;
//   const int K = 5;
//   const int S = pi.length();
//   const int Kpar = theta.length();
//
//
//   // Transition probs
//   Rcpp::NumericVector fn(T);
//   for(int t = 0; t < T; t++) {
//
//     Rcpp::NumericVector state0 = states_n( t, Rcpp::_ );
//     Rcpp::NumericVector state1 = states_n( t+1, Rcpp::_ );
//     int j = choices_n(t) - 1;
//
//     int k1 = state1( 3 ) - state0(3);
//     k1 = (k1 == 1000) ? 1 : (k1 == 2000) ? 2 : (k1 == 4000) ? 3 : 4 ;
//     int e1 = state1( 5 ) - 1;
//     int m1 = state1( 7 ) + 1; // this is (-1, 0, 1) -> (0, 1, 2)
//     int s1 = e1 * K * M  + m1*K + k1;
//
//     // // transition probabilities
//     Rcpp::NumericVector w_inc(E);
//     inc_transition_prob(w_inc, inc_transition_coefs, state0, E); // prob weights of income states
//
//     Rcpp::NumericMatrix w_savings(K, 2);
//     savings_transition_prob(w_savings, savings_transition_coefs, state0, j); // prob weights of savings states
//
//     Rcpp::NumericMatrix f_x(E*M*K,1);
//     match_transition_prob(
//       f_x,
//       match_transition_coefs,
//       state0,
//       j,
//       w_inc,
//       w_savings,
//       E,
//       M,
//       true
//     );
//     fn(t) = f_x(s1,0);
//
//   }
//
//
//   double L = 0;
//   Rcpp::NumericVector L_n(S);
//   Rcpp::NumericVector L_lambda_n(Kpar + S);
//   Rcpp::NumericMatrix H_L(Kpar+2*S,Kpar+2*S);
//   Rcpp::NumericMatrix grad_Lns_theta(S, Kpar) ; // L_n(s_n = s) gradients w.r.t. theta
//
//   Rcpp::NumericVector grad_Ln_p(S) ; // L_n gradient w.r.t. p(s)
//   Rcpp::NumericVector grad_p(S) ; // sum of ln L_t gradients w.r.t. p(s)
//
//   Rcpp::NumericMatrix H_lnL_theta_p(Kpar,S); // cross derivative of ln Lt w.r.t. theta and p(s)
//
//   // L_lambda_n(1:Kpar) = dlnL/dtheta
//   // exploit d[prod g(k)f(x,k)]/dx = d[exp(ln(prod g(k)f(x,k)))]/dx =
//   //   exp(ln(prod g(k)f(x,k))) * d[ln(prod g(k)f(x,k))]/dx =
//   //   prod g(k)f(x,k) * d[sum( ln(g(k)f(x,k)) )]/dx =
//   //   prod g(k)f(x,k) * sum( d[ln g(k) + ln f(x,k)]/dx ) =
//   //   prod g(k)f(x,k) [ sum( d[ln f(x,k)]/dx ) + sum(d[ln g(k)]/dx) ] =
//   //   prod g(k)f(x,k) sum( d[ln f(x,k)]/dx )
//
//   for(int s = 0; s < S; s++) {
//     L_n(s) = pi(s);
//     Rcpp::NumericVector grad_s(Kpar) ; // sum of ln L_t gradients
//     Rcpp::NumericMatrix H_theta_s(Kpar, Kpar) ; // sum of ln L_t hessians
//     // represents sum( d[ln f(x,k)]/dx )
//
//
//     for(int t = 0; t < T; t++) {
//       int index = t*J*S + s*J;
//       Rcpp::NumericMatrix v0 = v_n(
//         Rcpp::Range(index, index + J - 1), Rcpp::_
//       );
//       /*
//       Rcpp::NumericMatrix v0 = v_n(
//         Rcpp::Range(s*T*12 + t*12, s*T*12 + t*12 + 11), Rcpp::_
//       );
//       int p_diff_index = s*T*12 + t*12;
//        */
//
//       Rcpp::NumericVector prob(v0.nrow() + 1);
//       Choice_Prob(prob, v0, theta);
//
//
//       double Lt = prob(choices_n(t));
//       Lt *= fn(t);
//       L_n(s) *= Lt;
//
//       Gradient(grad_s, prob, choices_n(t), v0);
//       Hessian(H_theta_s, prob, v0);
//
//       Gradient_p(grad_p, prob, choices_n(t), pdiff_n, index, s);
//       Hessian_p_s(H_lnL_theta_p, prob, v0, pdiff_n, index, s);
//
//     }
//     L += L_n(s);
//
//     grad_Ln_p(s) = L_n(s) * grad_p(s);
//
//     for(int k = 0; k < Kpar; k++) {
//       grad_Lns_theta(s,k) = L_n(s) * grad_s(k); // Abl. der s-Komponente von L_lambda_n(k), * L
//       L_lambda_n(k) += grad_Lns_theta(s,k) ;
//       for(int l = k; l < Kpar; l++) {
//         H_L(k, l) += L_n(s) * (grad_s(k) * grad_s(l) + H_theta_s(k,l));
//       }
//       H_L(k, Kpar + s) += grad_Lns_theta(s,k) / pi(s); // Abl. von L_lambda_n(k) nach pi(s), * L
//
//       // w.r.t p(s)
//       H_L(k, Kpar + S + s) = L_n(s) * ( grad_s(k) * grad_p(s) + H_lnL_theta_p(k,s));
//
//     }
//   }
//
//
//   for(int k = 0; k < Kpar; k++) {
//     L_lambda_n(k) = L_lambda_n(k) / L;
//   }
//   for(int s = 0; s < S; s++) {
//     L_lambda_n(Kpar + s) = L_n(s) / (pi(s) * L);
//   }
//
//   for(int k = 0; k < Kpar; k++) {
//     for(int l = k; l < Kpar; l++) {
//       // L_lambda_n(theta(k)) w.r.t. theta(l)
//       H_L(k,l) = H_L(k,l) / L - L_lambda_n(k) * L_lambda_n(l);
//       // L_lambda_n(theta(l)) w.r.t. theta(k)
//       H_L(l,k) = H_L(k,l);
//     }
//     for(int s = 0; s < S; s++) {
//       // L_lambda_n(theta(k)) w.r.t. pi(s)
//       H_L(k, Kpar + s) = H_L(k, Kpar + s) / L - L_lambda_n(k) * L_lambda_n(Kpar+s) ;
//       // L_lambda_n(pi(s)) w.r.t. theta(k)
//       H_L(Kpar + s, k) = H_L(k, Kpar + s) ;  // Symmetrie wg. Satz von Schwarz
//
//       // Q_n w.r.t. theta
//       H_L(Kpar + S + s, k) = w_res_n(s) * (grad_Lns_theta(s,k)/L - q_n(s) * L_lambda_n(k)); // w_res_n(s) statt res_n(s) (!)
//
//       // L_lambda_n(theta) w.r.t. p
//       H_L(k, Kpar + S + s) = H_L(k, Kpar + S + s) / L - L_lambda_n(k) * grad_Ln_p(s) / L;
//
//     }
//   }
//
//   for(int s = 0; s < S; s++) {
//     for(int u = 0; u < S; u ++) {
//       // L_lambda_n(pi(s)) w.r.t. pi(u)
//       H_L(Kpar+s, Kpar+u) = -L_lambda_n(Kpar + s) * L_lambda_n(Kpar + u); // Abl. von L_lambda_n(Kpar+s) nach pi(u)
//
//       // Q_n(s) w.r.t. pi(u)
//       H_L(Kpar + S + s, Kpar + u) = -w_res_n(s) * q_n(u) * q_n(s) /pi(u) ;  // w_res_n(s) statt res_n(s) (!)
//
//       // L_lambda_n(pi(s)) w.r.t. p(u)
//       H_L(Kpar+s, Kpar + S + u) = -q_n(u) * q_n(s) /pi(u) * grad_Ln_p(s);
//
//       // Q_n(s) w.r.t. p(u)
//       H_L(Kpar + S + s, Kpar+S+u) = pi(s) * H_L(Kpar+s,Kpar+S+u) ;  // - off-diagonal entry
//     }
//
//     // Q_n(s) w.r.t. pi(s)
//     H_L(Kpar + S + s, Kpar + s) += w_res_n(s) * q_n(s) / pi(s)  ; // (!)
//
//     // L_lambda_n(pi(s)) w.r.t. p(s)
//     H_L(Kpar + s, Kpar + S + s) += q_n(s)/pi(s) * grad_Ln_p(s);
//
//     // Q_n(s) w.r.t. p(s)
//     H_L(Kpar + S + s, Kpar + S + s) += q_n(s) * w_n ; //  q_n(s) * w_n // (!) - diagonal entries
//   }
//
//   Rcpp::List L_out = Rcpp::List::create(L_lambda_n, H_L);
//   return L_out;
// }
//
//
//
//
//
// // /////////////// Simulation ///////////////
// //
// // [[Rcpp::export]]
// Rcpp::NumericVector update_state(
//   Rcpp::NumericVector& state,
//   const int& choice,
//   Rcpp::NumericVector& inc_transition_coefs,
//   Rcpp::NumericVector& inc_transition_errors,
//   Rcpp::NumericVector& sav_transition_coefs,
//   Rcpp::NumericVector& sav_transition_errors,
//   Rcpp::NumericMatrix& adu_transition_coefs,
//   Rcpp::NumericMatrix& kid_transition_coefs,
//   const double& p,
//   Rcpp::NumericVector& sizes_sqm
// ) {
//   const int j0 = choice - 1;
//
//   const int E = sav_transition_errors.length();
//
//   const int M = 2; // # of states adults variable
//   const int K = 3; // # of states kids variable
//
//   // transition probabilities (adults, kids)
//   Rcpp::NumericVector w_adu(M);
//   Rcpp::NumericVector w_kid(K);
//   // This does not depend on the choice.
//   adu_transition_prob(w_adu, adu_transition_coefs, state); // prob: adults in t+1
//   kid_transition_prob(w_kid, kid_transition_coefs, state);// prob: kids in t+1
//   Rcpp::NumericVector f_x (K * M); // prob weights (do not depend on the choice)
//
//   int m0 = M - 1;
//   int k0 = K - 1;
//   double p0 = 0.0;
//   for(int m = 0; m < M; m++) {
//     for(int k = 0; k < K; k++) {
//       p0 += w_adu(m) * w_kid(k);
//
//       if(p0 >= p) {
//         m0 = m;
//         k0 = k;
//         break;
//       }
//     }
//     if(p0 >= p) {
//       break ;
//     }
//   }
//
//   // income transition
//   Rcpp::NumericVector inc_lead (E);
//   inc_transition(
//     inc_lead,
//     inc_transition_coefs, inc_transition_errors,
//     state, m0
//   );
//
//   // savings transition (depends on income lead & choice)
//   Rcpp::NumericVector sav_lead (E);
//   sav_transition(
//     sav_lead,
//     sav_transition_coefs, sav_transition_errors,
//     state, inc_lead,
//     j0
//   );
//
//   Rcpp::NumericVector new_state( state.length() );
//
//   update_det_state(new_state, state, j0, sizes_sqm);
//   new_state(4) = sav_lead(0); //  imputed wealth
//   new_state(6) = inc_lead(0); //  income
//   new_state(7) = m0; // two_adults
//   new_state(8) = k0; // kids
//
//   new_state(19) = ((k0 > 0) & (state(8) > 0)) ? state(19) + 1 : 0 ; // age oldest child
//
//   return new_state;
// }
//
//
// //
// // # /*** R
// // # test()
// // # */


//
// // Alternative path: ######
// Rcpp::NumericMatrix mapped_dv_dr(
//     Rcpp::NumericVector& state,
//     const double& beta,
//     const int& J,
//     const int& terminal_j0,
//     const int& nvars_u_c,
//     const int& nvars_u_s,
//     const int& nvars_v,
//     const bool& compute_psi,
//     Rcpp::NumericVector& inc_transition_coefs,
//     Rcpp::NumericVector& inc_transition_errors,
//     Rcpp::NumericVector& sav_transition_coefs,
//     Rcpp::NumericVector& sav_transition_errors,
//     Rcpp::NumericMatrix& adu_transition_coefs,
//     Rcpp::NumericMatrix& kid_transition_coefs,
//     const int& hid,
//     const int& year,
//     const bool& save_W,
//     const bool& simulation,
//     const int& s0,
//     const int& S,
//     Rcpp::NumericVector& sizes_sqm,
//     Rcpp::NumericMatrix& rent_hs,
//     Rcpp::NumericVector& rent_hs_new
// ) {
//   const bool skip_hs = (rent_hs.ncol() == 1) ;
//
//   int nstate = state.length();
//
//   const int E = sav_transition_errors.length(); // error pairs from income and savings processes
//   const int M = 2; // # of states adults variable
//   const int K = 3; // # of states kids variable
//
//   const int nvars = nvars_u_c + nvars_u_s * S + nvars_v * 3;
//
//   const int terminal_j = terminal_j0 - 1;
//
//   Rcpp::NumericMatrix v (J * S, nvars + 2);
//
//   Rcpp::NumericMatrix States_j1 ((J-3) * E * M * K, nstate);
//   Rcpp::NumericMatrix States_01 (E * M * K, nstate);
//
//   // transition probabilities (adults, kids)
//   Rcpp::NumericVector w_adu(M);
//   Rcpp::NumericVector w_kid(K);
//   // This does not depend on the choice.
//   adu_transition_prob(w_adu, adu_transition_coefs, state); // prob: adults in t+1
//   kid_transition_prob(w_kid, kid_transition_coefs, state);// prob: kids in t+1
//   Rcpp::NumericVector f_x (E * M * K); // prob weights (do not depend on the choice)
//   for(int m = 0; m < M; m++) {
//     for(int k = 0; k < K; k++) {
//       for(int e = 0; e < E; e++) {
//         f_x(m * K * E + k * E + e) = w_adu(m) * w_kid(k) / E;
//       }
//     }
//   }
//
//   for(int j = 0; j < J; j++) {
//     // initialize containers
//     // Rcpp::NumericVector state_j1 (nstate); // deterministic state in t+1 if choice j
//     // Rcpp::NumericVector state_01 (nstate); // deterministic state in t+1 if choice 0
//     if( (j >= 0) & (j < 40)) {
//       du_dr(
//         dv_dr,
//         state, theta,
//         j, h,
//         sizes_sqm, rent_hs
//       );
//     }
//     mapped_vj(
//       v, // utility function vars (container)
//       States_j1, // states in t+1 if alternative choice in t (container -> psi)
//       States_01, // states in t+1 if baseline choice in t (container -> psi)
//       f_x, // container for prob weights
//       w_adu, // weights adults
//       w_kid, // weights kids
//       state, // current state
//       beta, // discount factor
//       j, // alternative j to be compared against baseline choice (j0 = -1)
//       terminal_j, // terminal choice in t+1
//       J, // total number of choices - 1
//       E, // number of potential income x savings states
//       M, // number of potential adult states
//       K, // number of potential kids states
//       S, // number of unobserved states
//       inc_transition_coefs, // coefficients for income transition function
//       inc_transition_errors, // errors income transition function
//       sav_transition_coefs, // coefficients for savings transition function
//       sav_transition_errors, // errors savings transition function
//       nvars,
//       nvars_u_c,
//       nvars_u_s,
//       nvars_v,
//       sizes_sqm,
//       rent_hs,
//       rent_hs_new,
//       skip_hs
//     );
//
//   }
//
//   Rcpp::NumericMatrix d_psi(J,40);
//   get_d_psi_sim(
//         d_psi,
//         States_j1,
//         States_01,
//         f_x,
//         S,
//         E,
//         K,
//         J,
//         M,
//         beta
//   );
//
//   return dv_dr;
// }
